Lichti et al., 2014 (Green Group) LC-MS Proteomics Re-analysis (Main effect of housing - EE vs IC rats NAc)

from ‘Environmental enrichment alters protein expression as well as the proteomic response to cocaine in rat nucleus accumbens’ in [frontiers in behavioral neuroscience] https://doi.org/10.3389/fnbeh.2014.00246

Also see Zhang et al., 2016 from Green Group https://doi.org/10.1016/j.neuroscience.2016.09.051

Contents

Load Packages

library(readxl) # read_excel()
library(tidyverse) # dplyr, stringr, ggplot2, readr, tidyr
library(gt) # for gt tables
library(UniProt.ws) # 1st test -  uniprot api 
library(org.Rn.eg.db) # Rat genome annotation package - incl. AnnotationDbi
library(clusterProfiler) #gsea
library(enrichplot) # for dotplot()

# library(pathview) # commented out atm | (but for visual kegg mapping)
# library(msigdbr) # commented out atm | MSigDB gene sets 
# # note some packages called directly instead e.g. cowplot:: & ggrepel::

Load initial Lichti df table

df_1 <- read_excel("input/Lichti_2014_datasheet_1.XLSX",
                   sheet = "quant-all",
                   trim_ws = TRUE) # does not trim white space behind accessions?
colnames(df_1)
##  [1] "Accession"                    "Name"                        
##  [3] "Gene Symbol"                  "EC log2 fold"                
##  [5] "EC est fold"                  "EC p value"                  
##  [7] "Coc log2 fold"                "Coc est fold"                
##  [9] "Coc p value"                  "interaction log2 effect size"
## [11] "interaction effect size"      "interaction p value"         
## [13] "fraction"
str(df_1)
## tibble [2,835 × 13] (S3: tbl_df/tbl/data.frame)
##  $ Accession                   : chr [1:2835] "Q3UHJ0" "Q3UHJ0" "P0C1X8" "P50475" ...
##  $ Name                        : chr [1:2835] "AAK1_MOUSE" "AAK1_MOUSE" "AAK1_RAT" "SYAC_RAT" ...
##  $ Gene Symbol                 : chr [1:2835] "AAK1" "AAK1" "Aak1" "AARS" ...
##  $ EC log2 fold                : num [1:2835] -0.442 0.404 0.216 -0.309 0.772 ...
##  $ EC est fold                 : num [1:2835] -1.36 1.32 1.16 -1.24 1.71 ...
##  $ EC p value                  : num [1:2835] 0.107 0.55 0.99 0.104 0.192 ...
##  $ Coc log2 fold               : num [1:2835] -0.312 -0.292 0.682 -0.15 -0.135 ...
##  $ Coc est fold                : num [1:2835] -1.24 -1.22 1.6 -1.11 -1.1 ...
##  $ Coc p value                 : num [1:2835] 0.226 0.679 0.99 0.396 0.885 ...
##  $ interaction log2 effect size: num [1:2835] 0.7 -0.357 0.195 0.442 -0.741 ...
##  $ interaction effect size     : num [1:2835] 1.62 -1.28 1.14 1.36 -1.67 ...
##  $ interaction p value         : num [1:2835] 0.0814 0.741 0.9899 0.1017 0.415 ...
##  $ fraction                    : chr [1:2835] "cytosolic" "nuclear" "membrane" "cytosolic" ...
length(unique(df_1$Accession)) # 1915
## [1] 1915
head(df_1, 10)
## # A tibble: 10 × 13
##    Accession Name        `Gene Symbol` `EC log2 fold` `EC est fold` `EC p value`
##    <chr>     <chr>       <chr>                  <dbl>         <dbl>        <dbl>
##  1 Q3UHJ0    AAK1_MOUSE  AAK1                 -0.442          -1.36       0.107 
##  2 Q3UHJ0    AAK1_MOUSE  AAK1                  0.404           1.32       0.550 
##  3 P0C1X8    AAK1_RAT    Aak1                  0.216           1.16       0.990 
##  4 P50475    SYAC_RAT    AARS                 -0.309          -1.24       0.104 
##  5 P50554    GABT_RAT    ABAT                  0.772           1.71       0.192 
##  6 P50554    GABT_RAT    ABAT                  0.0918          1.07       0.389 
##  7 P50554    GABT_RAT    ABAT                  0.670           1.59       0.990 
##  8 P06795    MDR1B_MOUSE Abcb1b               -2.35           -5.10       0.0906
##  9 Q5RKI8    ABCB8_RAT   ABCB8                -0.875          -1.83       0.0632
## 10 Q5RKI8    ABCB8_RAT   ABCB8                -0.0793         -1.06       0.991 
## # ℹ 7 more variables: `Coc log2 fold` <dbl>, `Coc est fold` <dbl>,
## #   `Coc p value` <dbl>, `interaction log2 effect size` <dbl>,
## #   `interaction effect size` <dbl>, `interaction p value` <dbl>,
## #   fraction <chr>
tail(df_1, 10)
## # A tibble: 10 × 13
##    Accession Name        `Gene Symbol` `EC log2 fold` `EC est fold` `EC p value`
##    <chr>     <chr>       <chr>                  <dbl>         <dbl>        <dbl>
##  1 P68511    1433F_RAT   YWHAH              -0.307            -1.24       0.0620
##  2 P68511    1433F_RAT   YWHAH               0.263             1.20       0.605 
##  3 P68511    1433F_RAT   YWHAH               0.000514          1.00       1.000 
##  4 P68254    1433T_MOUSE YWHAQ              -0.220            -1.17       0.255 
##  5 P68254    1433T_MOUSE YWHAQ              -0.296            -1.23       0.947 
##  6 P68254    1433T_MOUSE YWHAQ               0.491             1.41       0.990 
##  7 P63101    1433Z_MOUSE YWHAZ              -0.0697           -1.05       0.416 
##  8 P63101    1433Z_MOUSE YWHAZ               0.470             1.39       0.479 
##  9 P63101    1433Z_MOUSE YWHAZ               0.146             1.11       0.990 
## 10 Q8VIL3    ZWINT_RAT   Zwint              -0.131            -1.09       0.541 
## # ℹ 7 more variables: `Coc log2 fold` <dbl>, `Coc est fold` <dbl>,
## #   `Coc p value` <dbl>, `interaction log2 effect size` <dbl>,
## #   `interaction effect size` <dbl>, `interaction p value` <dbl>,
## #   fraction <chr>

Get extra uniprot in via uniprot.ws

Getting info like protein names, etc

ratUp <- UniProt.ws(10116) # rat
my_ids <- df_1$Accession
# https://www.uniprot.org/help/return_fields
# https://www.uniprot.org/help/general_annotation
# takes the rat & mouse uniprot accession IDs in Lichti's df & maps further uniprot info onto them
up_rat_mouse <- mapUniProt(
  from = "UniProtKB_AC-ID",
  to = "UniProtKB",
  columns = c("accession", "id", "reviewed", "protein_name", "gene_names",
              "gene_primary", "organism_name", "length", "keywordid", "keyword", "xref_geneid"),
  query = my_ids
  #query = list(taxId = 10116, ids = my_ids) # isn't selective to only rats?
)
# add website links
# follows a certain stucture (doesn't account all but has most)
up_rat_mouse_entry <- up_rat_mouse %>%
  mutate(uniprot_link = paste0("https://www.uniprot.org/uniprotkb/", Entry, "/entry")) %>% # add uniprot entry link! 
  mutate(wiki_link = paste0("https://en.wikipedia.org/wiki/", str_to_upper(Gene.Names..primary.))) # add wikipedia (needs 2 be uppercase)
  #most of the wiki ones work e.g. in if 4 or more characters
# write out basic df
# write_excel_csv(up_rat_mouse_entry, "output/2025_12_08_lichti_2014_test_map_v1.csv")

Bind uniprot info & df + add other info & clean up data

df_2 <- dplyr::left_join(df_1, up_rat_mouse_entry,
                         by = join_by(Accession == From)) %>%
  mutate(Direction = ifelse(`EC log2 fold` < 0, "Down", "Up"),
         gene_name = str_to_title(`Gene.Names..primary.`),
         .after = `Gene Symbol`) %>%
  mutate(Significant = ifelse(`EC p value` <= 0.05, 'yes', 'no'), .after = `EC p value`)  %>%
  filter(str_detect(Name, "RAT")) # 2835 to 1506
## Warning in dplyr::left_join(df_1, up_rat_mouse_entry, by = join_by(Accession == : Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 411 of `x` matches multiple rows in `y`.
## ℹ Row 1 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
##   "many-to-many"` to silence this warning.
df_3 <- df_2 %>%
  dplyr::select(Accession:Significant, fraction:wiki_link) %>% 
  mutate(across(where(is.numeric), ~ round(., 3))) %>% #if numeric col, round to 3 dp
  group_by(fraction, Direction) %>%  
  arrange(`EC p value`, `EC log2 fold`)

Plots

ggplot rat only non sig dot plot

# non significant plot - looks messy but you get the idea
ns_plot <- df_3 %>%
  ggplot(aes(x = gene_name, y = `EC log2 fold`, color = Direction, alpha = Significant)) +
  geom_point(shape = 1) +
  geom_hline(yintercept = 0, color = 'grey20') +
  scale_color_manual(breaks = c("Up", "Down"),
                     values = c("#dc322f", "#268bd2")) + 
  scale_alpha_manual(values = c(0.2,1)) +
  labs(title = 'Lichti et al., 2014 - EC log2fc x gene name - NAc',
       subtitle = 'Rats only (1506/2835), non-sig incl., split by fraction, LC-MS proteomics') +
  facet_wrap(~fraction, space = 'free_x', scale = 'free_x') + 
  theme_bw() 

ns_plot

ggplot lollipop chart (sig only)

# assign colors based on direction. | works by nvm ordering is weird
#label_colors <- ifelse(df_3$`EC log2 fold` > 0 , "#dc322f", "#268bd2")

# lollipop - ggplot ver - a big ugly
sig_plot <- df_3 %>%
  filter(`EC p value` <= 0.05) %>% # 139
  ggplot(aes(x = reorder(gene_name, `EC log2 fold`), # orders by log fold
                         y = `EC log2 fold`,
             color = Direction)) +
  geom_segment(aes(x = gene_name,
                   yend = `EC log2 fold`, y = 0)) +
  geom_point(aes(alpha = `EC p value`),
             shape = 19, size = 4) +
  scale_alpha(range = c(1,.1), guide = 'none') + # rev alpha so darker = more significant
  # geom_text(aes(label = gene_name), color = "grey10", size = 2.5,
  #           position = position_stack()) +
  geom_hline(yintercept = 0, color = 'grey20', linewidth = 0.1) +
  scale_color_manual(breaks = c("Up", "Down"),
                     values = c("#dc322f", "#268bd2"),
                     guide = 'none') + 
  labs(title = 'Lichti et al., 2014 - EC Log2FC x Gene Name - NAc',
       subtitle = 'Rats only (139/1506/2835), sig only, by fraction, LC-MS proteomics',
       y = 'EC Log2FC',
       x = 'Gene Name') +
  #coord_flip() + #fliptheplot
  facet_wrap(~fraction, space = 'free_x', scale = 'free_x') +
  #theme_void() +
  theme_bw() +
  theme(
    #axis.title = element_blank(),
    #axis.text.x = element_blank() 
    axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5, # 90, hj1 t/R, vj.5 cent
                               size = 9, face = 'italic'), 
    axis.text.y = element_text(hjust = 1, color = "grey20", size = 10), # rota 90 & hjust=1 (top/right)
        ) 

sig_plot

ggplot volcano

volcano_adjpv <- ggplot(df_3, 
                     aes(x =`EC log2 fold`,
                         y = -log2(`EC p value`), label = gene_name)) +
 geom_point(aes(color = case_when(
   `EC p value` < 0.05 & `EC log2 fold` > 0 ~ "Upregulated (p<0.05)",
   `EC p value` < 0.05 & `EC log2 fold` < 0 ~ "Downregulated (p<0.05)",
   `EC p value` < 0.1 & `EC log2 fold` > 0 ~ "Upregulated (p<0.1)",
   `EC p value` < 0.1 & `EC log2 fold` < 0 ~ "Downregulated (p<0.1)",
   TRUE ~ "Not Significant")),
   alpha = 1, size = 2.5, shape = 1) +
  xlim(-7.5,7.5) +
  geom_vline(xintercept = -1.5, linewidth =0.2, linetype = "dashed", color = "grey50") +
  geom_vline(xintercept = 1.5, linewidth =0.2, linetype = "dashed", color = "grey50") +
  geom_hline(yintercept = -log2(0.05), linetype = "dashed", linewidth = 0.2,  color = "grey30") +
  scale_color_manual(values = c("firebrick3", 
                                "royalblue4", 
                                "lightpink1", 
                                "lightblue3" , 
                                 "grey80"),
                     breaks = c("Upregulated (p<0.05)", 
                                "Downregulated (p<0.05)", 
                                "Upregulated (p<0.1)", 
                                "Downregulated (p<0.1)", 
                                "Not Significant"),
                     labels = c("Upregulated Proteins (up to p<0.05)", 
                                "Downregulated Proteins (up to p<0.05)", 
                                "Upregulated Proteins (up to p<0.1)", 
                                "Downregulated Proteins (up to p<0.1)" , 
                                "Not Significant")) +
  
  ggrepel::geom_text_repel() +
  facet_wrap(~ fraction, ncol = 1, scales = "free_y") +
  labs(
    title = "Lichti et al., 2014 LC-MS proteomics NAc rats - main effect of housing (EE vs IC) (2025_12_26)",
       subtitle = "Split by fraction | Volcano plots made from their provided df, not 100% sure on their p-value if adj or not?",
       x = "log2FC",
       y = "-Log2 (p-value)",
       color = "Protein Regulation") +
  geom_hline(yintercept = -log2(0.05), linetype = "dashed", linewidth = 0.2,  color = "grey30") +
  # annotate("text", x = .86, y = -log2(0.05), label = "Adj. p-value cutoff (0.05) \n",
  #          fontface = 'italic', size = 4, color = "grey10") +
  # geom_hline(yintercept = -log2(0.1), linetype = "dashed", linewidth = 0.2,  color = "grey50") +
  # annotate("text", x = .86, y = -log2(0.1), label = "Adj. p-value cutoff (0.1) \n",
  #          fontface = 'italic', size = 4, color = "grey30") +
  theme_bw() + 
  theme(legend.position = "bottom",
        strip.text = element_text(size = 12, face = "bold", color = "grey10")
        )


volcano_adjpv
## Warning: ggrepel: 626 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
## Warning: ggrepel: 331 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
## Warning: ggrepel: 441 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

gt Table

gt table for sig proteins only

df_3 |>
  dplyr::select(fraction, gene_name, `Protein.names`, Direction,
                `EC log2 fold`, `EC p value`, Entry) |> #cols for table
  dplyr::filter(`EC p value` <= 0.05) |>
  gt(
    groupname_col = dplyr::group_vars(df_3)
  ) |>
  tab_header(
    title = "Lichti et al., 2014 - EC vs. IC - NAc",
    subtitle = "Cleaned lc-ms proteomics dataframe, signficant rats only (139/1506/2835
    ), organised by fraction + direction & arranged by pvalue + EC log2FC"
  ) |>
   opt_stylize(style = 1, color = 'gray', add_row_striping = TRUE)
Lichti et al., 2014 - EC vs. IC - NAc
Cleaned lc-ms proteomics dataframe, signficant rats only (139/1506/2835 ), organised by fraction + direction & arranged by pvalue + EC log2FC
gene_name Protein.names EC log2 fold EC p value Entry
cytosolic - Up
Glod4 Glyoxalase domain-containing protein 4 0.538 0.007 Q5I0D1
Gsto1 Glutathione S-transferase omega-1 (GSTO-1) (EC 2.5.1.18) (Glutathione S-transferase omega 1-1) (GSTO 1-1) (Glutathione-dependent dehydroascorbate reductase) (EC 1.8.5.1) (Monomethylarsonic acid reductase) (MMA(V) reductase) (EC 1.20.4.2) (S-(Phenacyl)glutathione reductase) (SPG-R) 0.389 0.013 Q9Z339
Tpi1 Triosephosphate isomerase (TIM) (EC 5.3.1.1) (Methylglyoxal synthase) (EC 4.2.3.3) (Triose-phosphate isomerase) 0.274 0.016 P48500
Prdx5 Peroxiredoxin-5, mitochondrial (EC 1.11.1.24) (Antioxidant enzyme B166) (AOEB166) (PLP) (Peroxiredoxin V) (Prx-V) (Peroxisomal antioxidant enzyme) (Thioredoxin peroxidase PMP20) (Thioredoxin-dependent peroxiredoxin 5) 0.293 0.018 Q9R063
Arpc5l Actin-related protein 2/3 complex subunit 5-like protein (Arp2/3 complex 16 kDa subunit 2) (ARC16-2) 0.510 0.018 A1L108
Wdr7 WD repeat-containing protein 7 (TGF-beta resistance-associated protein TRAG) 0.652 0.018 Q9ERH3
Akr1a1 Aldo-keto reductase family 1 member A1 (EC 1.1.1.2) (EC 1.1.1.372) (EC 1.1.1.54) (3-DG-reducing enzyme) (Alcohol dehydrogenase [NADP(+)]) (Aldehyde reductase) (Glucuronate reductase) (EC 1.1.1.19) (Glucuronolactone reductase) (EC 1.1.1.20) (S-nitroso-CoA reductase) (ScorR) (EC 1.6.-.-) 0.258 0.020 P51635
Fahd2a Oxaloacetate tautomerase Fahd2a, mitochondrial (EC 5.3.2.2) (Fumarylacetoacetate hydrolase domain-containing protein 2A) 0.365 0.020 B2RYW9
Slc25a11 Mitochondrial 2-oxoglutarate/malate carrier protein (OGCP) (alpha-oxoglutarate carrier) (Solute carrier family 25 member 11) (SLC25A11) 0.430 0.020 P97700
Kyat3 Kynurenine--oxoglutarate transaminase 3 (EC 2.6.1.7) (Cysteine-S-conjugate beta-lyase 2) (EC 4.4.1.13) (Kynurenine aminotransferase 3) (Kynurenine aminotransferase III) (KATIII) (Kynurenine--glyoxylate transaminase) (EC 2.6.1.63) (Kynurenine--oxoglutarate transaminase III) 0.680 0.020 Q58FK9
Atp5f1b ATP synthase F(1) complex catalytic subunit beta, mitochondrial (EC 7.1.2.2) (ATP synthase F1 subunit beta) 0.161 0.021 P10719
Dnm1 Dynamin-1 (EC 3.6.5.5) (B-dynamin) (D100) (Dynamin I) (Dynamin, brain) 0.249 0.021 P21575
Dlat Dihydrolipoyllysine-residue acetyltransferase component of pyruvate dehydrogenase complex, mitochondrial (EC 2.3.1.12) (70 kDa mitochondrial autoantigen of primary biliary cirrhosis) (PBC) (Dihydrolipoamide acetyltransferase component of pyruvate dehydrogenase complex) (Pyruvate dehydrogenase complex component E2) (PDC-E2) (PDCE2) 0.273 0.021 P08461
Stip1 Stress-induced-phosphoprotein 1 (STI1) (Hsc70/Hsp90-organizing protein) (Hop) 0.299 0.021 O35814
Qdpr Dihydropteridine reductase (EC 1.5.1.34) (HDHPR) (Quinoid dihydropteridine reductase) 0.351 0.021 P11348
Septin3 Neuronal-specific septin-3 (G-septin) (P40) 0.571 0.021 Q9WU34
Slc4a10 Sodium-driven chloride bicarbonate exchanger (Solute carrier family 4 member 10) 3.018 0.021 Q80ZA5
Crmp1 Dihydropyrimidinase-related protein 1 (DRP-1) (Collapsin response mediator protein 1) (CRMP-1) (Inactive dihydropyrimidinase) 0.190 0.022 Q62950
Atp5po ATP synthase peripheral stalk subunit OSCP, mitochondrial (ATP synthase subunit O) (Oligomycin sensitivity conferral protein) (OSCP) (Sperm flagella protein 4) 0.394 0.022 Q06647
Nudt3 Diphosphoinositol polyphosphate phosphohydrolase 1 (DIPP-1) (EC 3.6.1.52) (Diadenosine hexaphosphate hydrolase) (Ap6A hydrolase) (EC 3.6.1.61) (Endopolyphosphatase) (EC 3.6.1.10) (Nucleoside diphosphate-linked moiety X motif 3) (Nudix motif 3) (m7GpppN-mRNA hydrolase) (EC 3.6.1.62) (m7GpppX diphosphatase) (EC 3.6.1.59) 0.793 0.025 Q566C7
Fahd1 Oxaloacetate tautomerase FAHD1, mitochondrial (EC 5.3.2.2) (Acylpyruvase FAHD1) (EC 3.7.1.5) (Fumarylacetoacetate hydrolase domain-containing protein 1) (Oxaloacetate decarboxylase) (OAA decarboxylase) (EC 4.1.1.112) 1.525 0.027 Q6AYQ8
Uba1 Ubiquitin-like modifier-activating enzyme 1 (EC 6.2.1.45) (Ubiquitin-activating enzyme E1) 0.197 0.028 Q5U300
Atp5f1c ATP synthase F(1) complex subunit gamma, mitochondrial (ATP synthase F1 subunit gamma) (F-ATPase gamma subunit) 0.384 0.028 P35435
Phgdh D-3-phosphoglycerate dehydrogenase (3-PGDH) (EC 1.1.1.95) 0.450 0.028 O08651
Mpst 3-mercaptopyruvate sulfurtransferase (MST) (EC 2.8.1.2) 0.332 0.029 P97532
Hibadh 3-hydroxyisobutyrate dehydrogenase, mitochondrial (HIBADH) (EC 1.1.1.31) 0.409 0.029 P29266
Nme1 Nucleoside diphosphate kinase A (NDK A) (NDP kinase A) (EC 2.7.4.6) (Metastasis inhibition factor NM23) (Tumor metastatic process-associated protein) 0.299 0.030 Q05982
Cryab Alpha-crystallin B chain (Alpha(B)-crystallin) 0.621 0.030 P23928
Hspa4 Heat shock 70 kDa protein 4 (Ischemia responsive 94 kDa protein) 0.200 0.031 O88600
Cs Citrate synthase, mitochondrial (EC 2.3.3.1) (Citrate (Si)-synthase) 0.322 0.031 Q8VHF5
Hspa8 Heat shock cognate 71 kDa protein (EC 3.6.4.10) (Heat shock 70 kDa protein 8) 0.208 0.032 P63018
Aco1 Cytoplasmic aconitate hydratase (Aconitase) (EC 4.2.1.3) (Citrate hydro-lyase) (Iron regulatory protein 1) (IRP1) (Iron-responsive element-binding protein 1) (IRE-BP 1) 0.271 0.032 Q63270
Atp5mf ATP synthase F(0) complex subunit f, mitochondrial (ATP synthase membrane subunit f) 0.332 0.032 D3ZAF6
Slc12a5 Solute carrier family 12 member 5 (Electroneutral potassium-chloride cotransporter 2) (Furosemide-sensitive K-Cl cotransporter) (K-Cl cotransporter 2) (rKCC2) (Neuronal K-Cl cotransporter) 0.370 0.032 Q63633
Slc6a17 Sodium-dependent neutral amino acid transporter SLC6A17 (Sodium-dependent neurotransmitter transporter NTT4) (Solute carrier family 6 member 17) 0.379 0.032 P31662
Dpp7 Dipeptidyl peptidase 2 (EC 3.4.14.2) (Dipeptidyl aminopeptidase II) (Dipeptidyl peptidase 7) (Dipeptidyl peptidase II) (DPP II) (Quiescent cell proline dipeptidase) 0.467 0.032 Q9EPB1
Atp5mg ATP synthase F(0) complex subunit g, mitochondrial (ATPase subunit g) (ATP synthase membrane subunit g) 0.599 0.032 Q6PDU7
Actrt1 Actin-related protein T1 1.752 0.032 Q5XIK1
Hsp90aa1 Heat shock protein HSP 90-alpha (EC 3.6.4.10) (Heat shock 86 kDa) (HSP 86) (HSP86) 0.336 0.033 P82995
Coq9 Ubiquinone biosynthesis protein COQ9, mitochondrial 1.045 0.033 Q68FT1
Dlg2 Disks large homolog 2 (Channel-associated protein of synapse-110) (Chapsyn-110) (Postsynaptic density protein PSD-93) 1.625 0.033 Q63622
Exoc2 Exocyst complex component 2 (Exocyst complex component Sec5) (rSec5) 1.769 0.033 O54921
Ppp5c Serine/threonine-protein phosphatase 5 (PP5) (EC 3.1.3.16) (Protein phosphatase T) (PPT) 2.036 0.033 P53042
Mpc2 Mitochondrial pyruvate carrier 2 (Brain protein 44) (Protein 0-44) 0.736 0.035 P38718
Ctsd Cathepsin D (EC 3.4.23.5) [Cleaved into: Cathepsin D 12 kDa light chain; Cathepsin D 9 kDa light chain; Cathepsin D 34 kDa heavy chain; Cathepsin D 30 kDa heavy chain] 1.318 0.035 P24268
Nfasc Neurofascin 0.172 0.036 P97685
Gnao1 Guanine nucleotide-binding protein G(o) subunit alpha (EC 3.6.5.-) 0.196 0.036 P59215
Pygb Glycogen phosphorylase, brain form (EC 2.4.1.1) 0.229 0.036 P53534
Tubb2a Tubulin beta-2A chain 0.316 0.036 P85108
Ppp2r2a Serine/threonine-protein phosphatase 2A 55 kDa regulatory subunit B alpha isoform (PP2A subunit B isoform B55-alpha) (B55) (PP2A subunit B isoform BRA) (PP2A subunit B isoform PR55-alpha) (PP2A subunit B isoform R2-alpha) (PP2A subunit B isoform alpha) 0.321 0.036 P36876
Ran GTP-binding nuclear protein Ran (EC 3.6.5.-) (GTPase Ran) (Ras-like protein TC4) (Ras-related nuclear protein) 0.505 0.036 P62828
Dpp10 Inactive dipeptidyl peptidase 10 (Dipeptidyl peptidase X) (DPP X) (Kv4 potassium channel auxiliary subunit) 0.579 0.036 Q6Q629
Vcp Transitional endoplasmic reticulum ATPase (TER ATPase) (EC 3.6.4.6) (15S Mg(2+)-ATPase p97 subunit) (Valosin-containing protein) (VCP) 0.184 0.037 P46462
Idh3b Isocitrate dehydrogenase [NAD] subunit beta, mitochondrial (Isocitric dehydrogenase subunit beta) (NAD(+)-specific ICDH subunit beta) 0.236 0.037 Q68FX0
Tubb3 Tubulin beta-3 chain (Neuron-specific class III beta-tubulin) 0.616 0.040 Q4QRB4
Usp7 Ubiquitin carboxyl-terminal hydrolase 7 (EC 3.4.19.12) (Deubiquitinating enzyme 7) (Herpesvirus-associated ubiquitin-specific protease) (rHAUSP) (Ubiquitin thioesterase 7) (Ubiquitin-specific-processing protease 7) 0.914 0.040 Q4VSI4
Mdh1 Malate dehydrogenase, cytoplasmic (EC 1.1.1.37) (Aromatic alpha-keto acid reductase) (KAR) (EC 1.1.1.96) (Cytosolic malate dehydrogenase) 0.279 0.041 O88989
Pgrmc1 Membrane-associated progesterone receptor component 1 (mPR) (25-DX) (Acidic 25 kDa protein) (Ventral midline antigen) (VEMA) 0.281 0.041 P70580
Fah Fumarylacetoacetase (FAA) (EC 3.7.1.2) (Beta-diketonase) (Fumarylacetoacetate hydrolase) 0.303 0.041 P25093
Idh1 Isocitrate dehydrogenase [NADP] cytoplasmic (IDH) (IDH1) (EC 1.1.1.42) (Cytosolic NADP-isocitrate dehydrogenase) (IDPc) (NADP(+)-specific ICDH) (Oxalosuccinate decarboxylase) 0.331 0.041 P41562
Cstb Cystatin-B (Cystatin-beta) (Liver thiol proteinase inhibitor) (Stefin-B) 0.444 0.041 P01041
Capza1 F-actin-capping protein subunit alpha-1 (CapZ alpha-1) 0.568 0.041 B2GUZ5
Slc3a2 Amino acid transporter heavy chain SLC3A2 (4F2 cell-surface antigen heavy chain) (4F2hc) (Solute carrier family 3 member 2) (CD antigen CD98) 0.188 0.042 Q794F9
Sncb Beta-synuclein (Phosphoneuroprotein 14) (PNP 14) 0.315 0.042 Q63754
Ola1 Obg-like ATPase 1 0.316 0.042 A0JPJ7
Lap3 Cytosol aminopeptidase (EC 3.4.11.1) (Cysteinylglycine-S-conjugate dipeptidase) (EC 3.4.13.23) (Leucine aminopeptidase 3) (LAP-3) (Leucyl aminopeptidase) (LAP) (Peptidase S) (Proline aminopeptidase) (EC 3.4.11.5) (Prolyl aminopeptidase) 0.474 0.042 Q68FS4
Inpp4a Inositol polyphosphate-4-phosphatase type I A (Inositol polyphosphate 4-phosphatase type I) (Type I inositol 3,4-bisphosphate 4-phosphatase) (EC 3.1.3.66) 0.225 0.043 Q62784
Pak1 Serine/threonine-protein kinase PAK 1 (EC 2.7.11.1) (Alpha-PAK) (Protein kinase MUK2) (p21-activated kinase 1) (PAK-1) (p68-PAK) 0.294 0.044 P35465
Gapdh Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) (EC 1.2.1.12) (38 kDa BFA-dependent ADP-ribosylation substrate) (BARS-38) (Peptidyl-cysteine S-nitrosylase GAPDH) (EC 2.6.99.-) 0.378 0.044 P04797
Cmpk1 UMP-CMP kinase (EC 2.7.4.14) (Deoxycytidylate kinase) (CK) (dCMP kinase) (Nucleoside-diphosphate kinase) (EC 2.7.4.6) (Uridine monophosphate/cytidine monophosphate kinase) (UMP/CMP kinase) (UMP/CMPK) 0.462 0.044 Q4KM73
Tkt Transketolase (TK) (EC 2.2.1.1) 0.192 0.045 P50137
Gna11 Guanine nucleotide-binding protein subunit alpha-11 (G alpha-11) (G-protein subunit alpha-11) 0.275 0.045 Q9JID2
Gpi Glucose-6-phosphate isomerase (GPI) (EC 5.3.1.9) (Autocrine motility factor) (AMF) (Neuroleukin) (NLK) (Phosphoglucose isomerase) (PGI) (Phosphohexose isomerase) (PHI) 0.161 0.046 Q6P6V0
Got1 Aspartate aminotransferase, cytoplasmic (cAspAT) (EC 2.6.1.1) (EC 2.6.1.3) (Cysteine aminotransferase, cytoplasmic) (Cysteine transaminase, cytoplasmic) (cCAT) (Glutamate oxaloacetate transaminase 1) (Transaminase A) 0.290 0.046 P13221
Psmb2 Proteasome subunit beta type-2 (Macropain subunit C7-I) (Multicatalytic endopeptidase complex subunit C7-I) (Proteasome component C7-I) (Proteasome subunit beta-4) (beta-4) 0.695 0.046 P40307
Tcp1 T-complex protein 1 subunit alpha (TCP-1-alpha) (EC 3.6.1.-) (CCT-alpha) 0.202 0.048 P28480
Tuba1b Tubulin alpha-1B chain (EC 3.6.5.-) (Alpha-tubulin 2) (Tubulin alpha-2 chain) [Cleaved into: Detyrosinated tubulin alpha-1B chain] 2.398 0.048 Q6P9V9
Sod2 Superoxide dismutase [Mn], mitochondrial (EC 1.15.1.1) 0.914 0.049 P07895
Akr1b1 Aldo-keto reductase family 1 member B1 (EC 1.1.1.21) (EC 1.1.1.300) (EC 1.1.1.372) (EC 1.1.1.54) (Aldehyde reductase) (Aldose reductase) (AR) 0.189 0.050 P07943
Pacsin1 Protein kinase C and casein kinase substrate in neurons protein 1 (Dynamin proline-rich domain-interacting protein) (Dynamin PRD-interacting protein) (Synaptic, dynamin-associated protein I) (Syndapin-1) (Syndapin-I) (SdpI) 0.189 0.050 Q9Z0W5
Prdx6 Peroxiredoxin-6 (EC 1.11.1.27) (1-Cys peroxiredoxin) (1-Cys PRX) (Acidic calcium-independent phospholipase A2) (aiPLA2) (EC 3.1.1.4) (Antioxidant protein 2) (Glutathione-dependent peroxiredoxin) (Lysophosphatidylcholine acyltransferase 5) (LPC acyltransferase 5) (LPCAT-5) (Lyso-PC acyltransferase 5) (EC 2.3.1.23) (Non-selenium glutathione peroxidase) (NSGPx) (Thiol-specific antioxidant protein) 0.195 0.050 O35244
Sfxn3 Sideroflexin-3 0.272 0.050 Q9JHY2
Bdh1 D-beta-hydroxybutyrate dehydrogenase, mitochondrial (EC 1.1.1.30) (3-hydroxybutyrate dehydrogenase) (BDH) 0.343 0.050 P29147
Hnrnpd Heterogeneous nuclear ribonucleoprotein D0 (hnRNP D0) (AU-rich element RNA-binding protein 1) 0.856 0.050 Q9JJ54
cytosolic - Down
Arf4 ADP-ribosylation factor 4 -2.212 0.013 P61751
Adgrl1 Adhesion G protein-coupled receptor L1 (Calcium-independent alpha-latrotoxin receptor 1) (CIRL-1) (Latrophilin-1) -3.426 0.014 O88917
C1qbp Complement component 1 Q subcomponent-binding protein, mitochondrial (GC1q-R protein) (Glycoprotein gC1qBP) (C1qBP) -2.686 0.014 O35796
Cygb Cytoglobin (Nitric oxygen dioxygenase CYGB) (EC 1.14.12.-) (Nitrite reductase CYGB) (EC 1.7.-.-) (Pseudoperoxidase CYGB) (EC 1.11.1.-) (Stellate cell activation-associated protein) (Superoxide dismutase CYGB) (EC 1.15.1.1) -3.009 0.018 Q921A4
Gyg1 Glycogenin-1 (GN-1) (GN1) (EC 2.4.1.186) -2.965 0.018 O08730
Entpd2 Ectonucleoside triphosphate diphosphohydrolase 2 (NTPDase 2) (EC 3.6.1.-) (CD39 antigen-like 1) (Ecto-ATP diphosphohydrolase 2) (Ecto-ATPDase 2) (Ecto-ATPase 2) -2.420 0.018 O35795
Map2k2 Dual specificity mitogen-activated protein kinase kinase 2 (MAP kinase kinase 2) (MAPKK 2) (EC 2.7.12.2) (ERK activator kinase 2) (MAPK/ERK kinase 2) (MEK 2) -2.246 0.018 P36506
Vps35l VPS35 endosomal protein-sorting factor-like -0.910 0.018 Q5XI83
Csad Cysteine sulfinic acid decarboxylase (EC 4.1.1.29) (Aspartate 1-decarboxylase) (EC 4.1.1.11) (Cysteine-sulfinate decarboxylase) (Sulfinoalanine decarboxylase) -5.304 0.020 Q64611
Psmb7 Proteasome subunit beta type-7 (EC 3.4.25.1) (Macropain chain Z) (Multicatalytic endopeptidase complex chain Z) (Proteasome subunit Z) (Proteasome subunit beta-2) (beta-2) -3.307 0.020 Q9JHW0
Pnpo Pyridoxine-5'-phosphate oxidase (EC 1.4.3.5) (Pyridoxamine-phosphate oxidase) -3.138 0.020 O88794
Rabggta Geranylgeranyl transferase type-2 subunit alpha (EC 2.5.1.60) (Geranylgeranyl transferase type II subunit alpha) (Rab geranyl-geranyltransferase subunit alpha) (Rab GG transferase alpha) (Rab GGTase alpha) (Rab geranylgeranyltransferase subunit alpha) -1.157 0.020 Q08602
Ddx1 ATP-dependent RNA helicase DDX1 (EC 3.6.4.13) (DEAD box protein 1) -0.622 0.020 Q641Y8
Copb1 Coatomer subunit beta (Beta-coat protein) (Beta-COP) -0.455 0.020 P23514
Stx1a Syntaxin-1A (Neuron-specific antigen HPC-1) (Synaptotagmin-associated 35 kDa protein) (P35A) -0.431 0.025 P32851
Serpina1 Alpha-1-antiproteinase (Alpha-1-antitrypsin) (Alpha-1-proteinase inhibitor) (Serpin A1) -2.763 0.028 P17475
Fkbp1a Peptidyl-prolyl cis-trans isomerase FKBP1A (PPIase FKBP1A) (EC 5.2.1.8) (12 kDa FK506-binding protein) (12 kDa FKBP) (FKBP-12) (FK506-binding protein 1A) (FKBP-1A) (Immunophilin FKBP12) (Rotamase) -0.952 0.030 Q62658
Adk Adenosine kinase (AK) (EC 2.7.1.20) (Adenosine 5'-phosphotransferase) -2.509 0.031 Q64640
Man2c1 Alpha-mannosidase 2C1 (EC 3.2.1.24) (Alpha-D-mannoside mannohydrolase) (Mannosidase alpha class 2C member 1) -2.007 0.031 P21139
Por NADPH--cytochrome P450 reductase (CPR) (P450R) (EC 1.6.2.4) -1.797 0.031 P00388
Lcmt1 Leucine carboxyl methyltransferase 1 (EC 2.1.1.233) ([Phosphatase 2A protein]-leucine-carboxy methyltransferase 1) -1.726 0.031 Q6P4Z6
Ppm1e Protein phosphatase 1E (EC 3.1.3.16) (Ca(2+)/calmodulin-dependent protein kinase phosphatase N) (CaMKP-N) (CaMKP-nucleus) (CaMKN) (Partner of PIX 1) (Partner of PIX-alpha) (Partner of PIXA) -1.311 0.031 Q80Z30
Pafah1b2 Platelet-activating factor acetylhydrolase IB subunit alpha2 (EC 3.1.1.47) (PAF acetylhydrolase 30 kDa subunit) (PAF-AH 30 kDa subunit) (PAF-AH subunit beta) (PAFAH subunit beta) (Platelet-activating factor acetylhydrolase alpha 2 subunit) (PAF-AH alpha 2) -0.666 0.031 O35264
Rpn2 Dolichyl-diphosphooligosaccharide--protein glycosyltransferase subunit 2 (Dolichyl-diphosphooligosaccharide--protein glycosyltransferase 63 kDa subunit) (Ribophorin II) (RPN-II) (Ribophorin-2) -0.681 0.033 P25235
Rala Ras-related protein Ral-A (EC 3.6.5.2) -1.185 0.035 P63322
Prepl Prolyl endopeptidase-like (EC 3.4.21.-) (Prolylendopeptidase-like) -0.682 0.035 Q5HZA6
Plcb1 1-phosphatidylinositol 4,5-bisphosphate phosphodiesterase beta-1 (EC 3.1.4.11) (PLC-154) (Phosphoinositide phospholipase C-beta-1) (Phospholipase C-I) (PLC-I) (Phospholipase C-beta-1) (PLC-beta-1) -0.253 0.039 P10687
Prune1 Exopolyphosphatase PRUNE1 (EC 3.6.1.1) -0.557 0.041 Q6AYG3
Slc8a1 Sodium/calcium exchanger 1 (Na(+)/Ca(2+)-exchange protein 1) (Solute carrier family 8 member 1) -0.490 0.041 Q01728
Nono Non-POU domain-containing octamer-binding protein (NonO protein) -2.858 0.043 Q5FVM4
Pir Pirin (EC 1.13.11.24) (Probable quercetin 2,3-dioxygenase PIR) (Probable quercetinase) -0.444 0.045 Q5M827
Sacm1l Phosphatidylinositol-3-phosphatase SAC1 (EC 3.1.3.64) (Phosphatidylinositol-4-phosphate phosphatase) (Suppressor of actin mutations 1-like protein) -1.268 0.050 Q9ES21
Ndufv2 NADH dehydrogenase [ubiquinone] flavoprotein 2, mitochondrial (EC 7.1.1.2) (NADH-ubiquinone oxidoreductase 24 kDa subunit) -0.608 0.050 P19234
nuclear - Up
Atp5me ATP synthase F(0) complex subunit e, mitochondrial (ATPase subunit e) (ATP synthase membrane subunit e) 1.618 0.017 P29419
Phb1 Prohibitin 1 1.098 0.022 P67779
Uqcrc1 Cytochrome b-c1 complex subunit 1, mitochondrial (Complex III subunit 1) (Core protein I) (Ubiquinol-cytochrome-c reductase complex core protein 1) 0.853 0.023 Q68FY0
Atp5po ATP synthase peripheral stalk subunit OSCP, mitochondrial (ATP synthase subunit O) (Oligomycin sensitivity conferral protein) (OSCP) (Sperm flagella protein 4) 0.953 0.023 Q06647
Timm13 Mitochondrial import inner membrane translocase subunit Tim13 0.708 0.034 P62076
Ndufa9 NADH dehydrogenase [ubiquinone] 1 alpha subcomplex subunit 9, mitochondrial (Complex I-39kD) (CI-39kD) (NADH-ubiquinone oxidoreductase 39 kDa subunit) (Sperm flagella protein 3) 0.561 0.045 Q5BK63
Atp5f1a ATP synthase F(1) complex subunit alpha, mitochondrial (ATP synthase F1 subunit alpha) 0.748 0.045 P15999
Phb2 Prohibitin-2 (B-cell receptor-associated protein BAP37) (BAP-37) 0.751 0.045 Q5XIH7
Uqcrc2 Cytochrome b-c1 complex subunit 2, mitochondrial (Complex III subunit 2) (Core protein II) (Ubiquinol-cytochrome-c reductase complex core protein 2) 0.909 0.045 P32551
Ndufv2 NADH dehydrogenase [ubiquinone] flavoprotein 2, mitochondrial (EC 7.1.1.2) (NADH-ubiquinone oxidoreductase 24 kDa subunit) 0.948 0.045 P19234
Uqcrfs1 Cytochrome b-c1 complex subunit Rieske, mitochondrial (EC 7.1.1.8) (Complex III subunit 5) (Cytochrome b-c1 complex subunit 5) (Liver regeneration-related protein LRRGT00195) (Rieske iron-sulfur protein) (RISP) (Rieske protein UQCRFS1) (Ubiquinol-cytochrome c reductase iron-sulfur subunit) [Cleaved into: Cytochrome b-c1 complex subunit 9 (Su9) (Subunit 9) (8 kDa subunit 9) (Complex III subunit IX) (Cytochrome b-c1 complex subunit 11) (UQCRFS1 mitochondrial targeting sequence) (UQCRFS1 MTS) (Ubiquinol-cytochrome c reductase 8 kDa protein)] 0.989 0.045 P20788
Ppia Peptidyl-prolyl cis-trans isomerase A (PPIase A) (EC 5.2.1.8) (Cyclophilin A) (Cyclosporin A-binding protein) (Rotamase A) (p1B15) (p31) [Cleaved into: Peptidyl-prolyl cis-trans isomerase A, N-terminally processed] 1.167 0.045 P10111
Cd47 Leukocyte surface antigen CD47 (Integrin-associated protein) (IAP) (CD antigen CD47) 0.538 0.047 P97829
Stxbp1 Syntaxin-binding protein 1 (N-Sec1) (Protein unc-18 homolog 1) (Unc18-1) (Protein unc-18 homolog A) (Unc-18A) (p67) (rbSec1) 0.663 0.047 P61765
Aldoa Fructose-bisphosphate aldolase A (EC 4.1.2.13) (Muscle-type aldolase) 0.675 0.047 P05065
Tomm22 Mitochondrial import receptor subunit TOM22 homolog (rTOM22) (Translocase of outer membrane 22 kDa subunit homolog) 0.796 0.047 Q75Q41
Maoa Amine oxidase [flavin-containing] A (EC 1.4.3.21) (EC 1.4.3.4) (Monoamine oxidase type A) (MAO-A) 0.818 0.047 P21396
Cox5a Cytochrome c oxidase subunit 5A, mitochondrial (Cytochrome c oxidase polypeptide Va) 0.833 0.047 P11240
Cs Citrate synthase, mitochondrial (EC 2.3.3.1) (Citrate (Si)-synthase) 0.900 0.047 Q8VHF5
Pkm Pyruvate kinase PKM (EC 2.7.1.40) (Pyruvate kinase muscle isozyme) (Threonine-protein kinase PKM2) (EC 2.7.11.1) (Tyrosine-protein kinase PKM2) (EC 2.7.10.2) 0.908 0.047 P11980
Atp5f1b ATP synthase F(1) complex catalytic subunit beta, mitochondrial (EC 7.1.2.2) (ATP synthase F1 subunit beta) 0.930 0.047 P10719
Aldoc Fructose-bisphosphate aldolase C (EC 4.1.2.13) (Brain-type aldolase) 1.140 0.047 P09117
Hbb Hemoglobin subunit beta-1 (Beta-1-globin) (Hemoglobin beta chain, major-form) (Hemoglobin beta-1 chain) 1.286 0.047 P02091
Cisd1 CDGSH iron-sulfur domain-containing protein 1 (Cysteine transaminase CISD1) (EC 2.6.1.3) (MitoNEET) 1.344 0.047 B0K020
Arhgdia Rho GDP-dissociation inhibitor 1 (Rho GDI 1) (Rho-GDI alpha) 1.765 0.047 Q5XI73
Hsph1 Heat shock protein 105 kDa (Heat shock 110 kDa protein) 0.427 0.050 Q66HA8
nuclear - Down
Dlg4 Disks large homolog 4 (Postsynaptic density protein 95) (PSD-95) (Synapse-associated protein 90) (SAP-90) (SAP90) -0.378 0.047 P31016
Shank3 SH3 and multiple ankyrin repeat domains protein 3 (Shank3) (Proline-rich synapse-associated protein 2) (ProSAP2) (SPANK-2) -0.439 0.049 Q9JLU4

Bioinformatics

old nooby 1 by 1 method (commented out)

Ranked List

already partially prepped from above

# # ranked combined list
# geneList_df <- df_3 %>%
#   mutate(ranking = sign(`EC log2 fold`)*(-log10(`EC p value`))) %>%
#   group_by(Entry) %>%
#   summarise(ranking = mean(ranking, na.rm = TRUE)) %>%
#   arrange(-ranking)
# 
# # convert 2 a named vector to supply to fgsea()
# geneList <- tibble::deframe(geneList_df) # 955
# 
# plot(geneList)
# # ranked cytosolic list
# geneList_df_c <- df_3 %>%
#   filter(fraction == 'cytosolic') %>%
#   mutate(ranking = sign(`EC log2 fold`)*(-log10(`EC p value`))) %>%
#   group_by(Entry) %>%
#   summarise(ranking = mean(ranking, na.rm = TRUE)) %>%
#   arrange(-ranking)
# 
# # convert 2 a named vector to supply to fgsea()
# geneList_c <- tibble::deframe(geneList_df_c) # 689
# 
# plot(geneList_c)
# # ranked membrane list
# geneList_df_m <- df_3 %>%
#   filter(fraction == 'membrane') %>%
#   mutate(ranking = sign(`EC log2 fold`)*(-log10(`EC p value`))) %>%
#   group_by(Entry) %>%
#   summarise(ranking = mean(ranking, na.rm = TRUE)) %>%
#   arrange(-ranking)
# 
# # convert 2 a named vector to supply to fgsea()
# geneList_m <- tibble::deframe(geneList_df_m) # 340
# 
# plot(geneList_m)
# # ranked nuclear list
# geneList_df_n <- df_3 %>%
#   filter(fraction == 'nuclear') %>%
#   mutate(ranking = sign(`EC log2 fold`)*(-log10(`EC p value`))) %>%
#   group_by(Entry) %>%
#   summarise(ranking = mean(ranking, na.rm = TRUE)) %>%
#   arrange(-ranking)
# 
# # convert 2 a named vector to supply to fgsea()
# geneList_n <- tibble::deframe(geneList_df_n) # 477
# 
# plot(geneList_n)

KEGG

Nuclear

# # nuclear
# kegg_res_n <- gseKEGG(geneList = geneList_n,
#                     organism = 'rno',
#                     keyType = 'uniprot',
#                     by = 'fgsea',
#                     pvalueCutoff = 0.05,
#                     pAdjustMethod = 'BH', #BejaminiHoc Correciton
#                     )
# 
# gse_kegg_tib_n <- as_tibble(kegg_res_n) # 21 sig
# enrichplot::dotplot(kegg_res_n,
#                     split = ".sign",  # split by + vs - enrirchment score
#                     showCategory = 30,
#                     title = "Nuclear - NAc - Lichti et al., 2014 (EE vs IC) - gseKEGG_r - adj.p<0.05") +
#   facet_grid(.~.sign) # split plot by seign)

Cytosolic

# # cytosolic
# kegg_res_c <- gseKEGG(geneList = geneList_c,
#                     organism = 'rno',
#                     keyType = 'uniprot',
#                     by = 'fgsea',
#                     pvalueCutoff = 0.05,
#                     pAdjustMethod = 'BH', #BejaminiHoc Correciton
#                     )
# 
# gse_kegg_tib_c <- as_tibble(kegg_res_c) # 21 sig
# enrichplot::dotplot(kegg_res_c,
#                     split = ".sign",  # split by + vs - enrirchment score
#                     showCategory = 30,
#                     title = "Cytosolic - NAc - Lichti et al., 2014 (EE vs IC) - gseKEGG_r - adj.p<0.05") +
#   facet_grid(.~.sign) # split plot by seign)

Membrane

# # membrane
# kegg_res_m <- gseKEGG(geneList = geneList_m,
#                     organism = 'rno',
#                     keyType = 'uniprot',
#                     by = 'fgsea',
#                     pvalueCutoff = 0.05,
#                     pAdjustMethod = 'BH', #BejaminiHoc Correciton
#                     )
# 
# gse_kegg_tib_m <- as_tibble(kegg_res_m) # 3 sig
# enrichplot::dotplot(kegg_res_m,
#                     split = ".sign",  # split by + vs - enrirchment score
#                     showCategory = 30,
#                     title = "Membrane - NAc - Lichti et al., 2014 (EE vs IC) - gseKEGG_r - adj.p<0.05") +
#   facet_grid(.~.sign) # split plot by seign)
# # oxidative phosphorylation - rno00190
# rno00190 <- pathview(gene.data = geneList_n,
#                      gene.idtype = "SYMBOL",
#                      pathway.id = 'rno00190', # oxidative phosphorylation - Rat
#                      species = 'rno',
#                      limit = list(gene=max(abs(geneList)), cpd=1))

Trying for loops for gsea (KEGG)

Setup containers & prepare info for sequences

# make vector of fractions (cytosolic, membrane, & nuclear) - in prep of sequences
fraction_names <- unique(df_3$fraction)

# make container lists to store results
# gsea_results_list <- vector(mode = "list", length = 10)
gsea_results_list <- vector(mode = "list") # actually make empty

Operation

for (fr in fraction_names){ # SEQUENCE - for each item (fr) in fraction_names
  
  # OPERATIONS
  # Part 1
  # make a ranked gene list for this fraciton
  geneList_working_df <- df_3 %>%
      # filter by item | fraction is column in df_3, fr is i (current item in seq)
    filter(fraction == fr) %>% 
      # compute ranking metric
    mutate(ranking = sign(`EC log2 fold`)*(-log10(`EC p value`))) %>%
      # for any duplicate entries, # average & remove any nas
    group_by(Entry) %>% 
    summarise(ranking = mean(ranking, na.rm = TRUE)) %>%
      # lastly, arrange from pos to neg 
    arrange(-ranking)
  
  # store in list (not necessary...)
  gsea_results_list$geneList_df[[fr]] <- geneList_working_df
  
  # Part 1.2
  # next, deframe into a named vector needed to supply to fgsea()
  geneList_vec_working <- tibble::deframe(geneList_working_df)
  
  # store in list 
  gsea_results_list$geneList_vec[[fr]] <- geneList_vec_working
  
  # Part 1.3 | plots fine but does not save?!?!?!
  # Plot ranked gene list
  gL_plot <- plot(geneList_vec_working, main = paste(fr, '- Ranked gene list'))
  print(gL_plot)

  gsea_results_list$geneList_plots[[fr]] <- gL_plot
  
 
  # Part 2
  # run KEGG gsea with clusterProfiler
  kegg_res_working <- clusterProfiler::gseKEGG(geneList = geneList_vec_working,
                              organism = 'rno',
                              keyType = 'uniprot',
                              by = 'fgsea',
                              pvalueCutoff = 0.05,
                              pAdjustMethod = 'BH', #BejaminiHoc Correciton
                              )
  
  # save kegg result original to list vec
  gsea_results_list$kegg_results[[fr]] <- kegg_res_working
  
  # not needed here so compute & save directly? or can do anyway? idk?
  gsea_results_list$kegg_results_tib[[fr]] <- tibble::as_tibble(kegg_res_working) 
  
  
  # Part 3
  # Make dotplots | this one plots & saves properly!
  # 3.1 - regular dot blot
  dp <- enrichplot::dotplot(kegg_res_working,
                            showCategory = 30,
                            font.size = 11,
                            title = paste0(stringr::str_glue(
                              '{fr} - NAc - Lichti, 2014 (EE vs IC) - gseKEGG_r'))) 

  
  print(dp)   # print dot plots 
  gsea_results_list$dotplots[[fr]] <- dp # save to list 
  
  # 3.2 - dot plot split by sign
  dp_sign <- enrichplot::dotplot(kegg_res_working,
                            split = '.sign',
                            showCategory = 30,
                            font.size = 11,
                            title = paste0(stringr::str_glue(
                              '{fr} - NAc - Lichti, 2014 (EE vs IC) - gseKEGG_r split by sign'))) +
    facet_grid(.~.sign)
  
  print(dp_sign) # print dot plots 
  gsea_results_list$dotplots_sign[[fr]] <- dp_sign # save to list

    
}

## NULL
## Reading KEGG annotation online: "https://rest.kegg.jp/link/rno/pathway"...
## Reading KEGG annotation online: "https://rest.kegg.jp/list/pathway/rno"...
## Reading KEGG annotation online: "https://rest.kegg.jp/conv/uniprot/rno"...
## using 'fgsea' for GSEA analysis, please cite Korotkevich et al (2019).
## preparing geneSet collections...
## GSEA analysis...
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (33.82% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
## leading edge analysis...
## done...

## NULL
## using 'fgsea' for GSEA analysis, please cite Korotkevich et al (2019).
## preparing geneSet collections...
## GSEA analysis...
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (30.61% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize =
## minSize, : For some of the pathways the P-values were likely overestimated. For
## such pathways log2err is set to NA.
## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize =
## minSize, : For some pathways, in reality P-values are less than 1e-10. You can
## set the `eps` argument to zero for better estimation.
## leading edge analysis...
## done...

## NULL
## using 'fgsea' for GSEA analysis, please cite Korotkevich et al (2019).
## preparing geneSet collections...
## GSEA analysis...
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (89.71% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
## leading edge analysis...
## done...

# horrific name but at least clear...
#Lichti_2014_gsea_results_2025_12_12_v1 <- gsea_results_list 

# save(Lichti_2014_gsea_results_2025_12_12_v1, 
#      file = "output/2025_12_12_Lichti_2014_gsea_KEGG_results_rats_NAc_v1.rda")
# Warning means that some genes have tied rankings
# note the membrane one is especially high 89.71%  - this corresponds to the weird # pattern in geneList ranking plot -  

# Warning messages:
# 1: In preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam,  :
#   There are ties in the preranked stats (33.82% of the list).
# The order of those tied genes will be arbitrary, which may produce unexpected results.
# 2: In preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam,  :
#   There are ties in the preranked stats (30.61% of the list).
# The order of those tied genes will be arbitrary, which may produce unexpected results.
# 3: In fgseaMultilevel(pathways = pathways, stats = stats, minSize = minSize,  :
#   For some of the pathways the P-values were likely overestimated. For such pathways log2err is set to NA.
# 4: In fgseaMultilevel(pathways = pathways, stats = stats, minSize = minSize,  :
#   For some pathways, in reality P-values are less than 1e-10. You can set the `eps` argument to zero for better estimation.
# 5: In preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam,  :
#   There are ties in the preranked stats (89.71% of the list).
# The order of those tied genes will be arbitrary, which may produce unexpected results.
check_1 <- gsea_results_list$geneList_df$cytosolic

length(unique(gsea_results_list$geneList_df$cytosolic$Entry)) # 689 unique proteins
## [1] 689
length(unique(gsea_results_list$geneList_df$cytosolic$ranking)) # 456 unique ranks
## [1] 456
# cow plot for combined + plotlist = to plot a list of plots 
cowplot::plot_grid(plotlist = c(gsea_results_list$dotplots,
                                gsea_results_list$dotplots_sign),
                   ncol = 2)

# print(gsea_results_list$dotplots)
# weird... highlights importance of plotting data!!!
plot(df_3$`EC log2 fold`)

2025_12_12 - Note to self maybe add some testing & progress of loop info - ref to epi cook book. Also try nested for loops with GO analysis e.g. for each of the 3 fractions (cytosolic, nuclear, membrane) run each of the 3 GO (CC, BP, MF)

Literature Review Table

Make lit review table (similar to McNair one)

colnames(df_3) # check col names 
##  [1] "Accession"            "Name"                 "Gene Symbol"         
##  [4] "Direction"            "gene_name"            "EC log2 fold"        
##  [7] "EC est fold"          "EC p value"           "Significant"         
## [10] "fraction"             "Entry"                "Entry.Name"          
## [13] "Reviewed"             "Protein.names"        "Gene.Names"          
## [16] "Gene.Names..primary." "Organism"             "Length"              
## [19] "Keyword.ID"           "Keywords"             "GeneID"              
## [22] "uniprot_link"         "wiki_link"
# lit review table 
lit_re <- df_3 %>%
  mutate(
    Space = " ", # empty space column
    Year = 2014,
    Author = 'Lichti',
    Animal = 'Rat',
    Condition = paste0("Fraction - ", fraction),
    Region = 'Nucleus_Accumbens',
    Whole = 'No',
    Exp_type = 'Protein',
    Method = 'PROTEOMICS_LCMS_LFQ_OS',
    Used_name = `Gene Symbol`,
    Mod = " ",
    Mod_type = " ",
    Gene_name = gene_name, # lazy way to avoid moving rows (don't need headers anyway)
    Comparison = 1,
    In_EE = Direction, # note was coded as >0 Up, <0 Down. Maybe change one to same?? 
    Fold_change = `EC log2 fold`,
    #Arrows = ifelse(Direction == 'Up', '+', '-'),
    Arrows = case_when( # dplyr case when, instead of ifelse
      `EC log2 fold` > 2.5 ~ '+++',
      `EC log2 fold` > 1 ~ '++',
      `EC log2 fold` > .05 ~ '+',
      `EC log2 fold` > -.05 ~ '=',
      `EC log2 fold` > -1 ~ '-',
      `EC log2 fold` > -2.5 ~ '--',
      `EC log2 fold` <= -2.5 ~ '---',
      ),
    #Significant = ifelse(`EC p value` <= 0.05, 'yes', 'no'),
    signficant = Significant,
    pvalue = `EC p value`,
    N_tested = " ",
    Stats_adj = 'adjusted',
    Note = ' ',
    Full_name = Protein.names,
    Uniprot_rat_accession = Entry,
    Uniprot_entry_name = Entry.Name,
    Uniprot_link = uniprot_link,
    Wikipedia_link = wiki_link,
    Curation_method = ifelse(!is.na(Entry), 'auto - uniprot_acc_query', 'manual'),
    Class_group = ' ',
    Other_link = ' ',
    Other_comments = '*misc col shows the col names in their excel df',
    misc = paste0(Accession, ";", `Gene Symbol`, ";", Name)
  )

Export

#write out | na values = blank 
# write_excel_csv(lit_re,
#                 "output/2025_12_12_Lichti_2014_rats_NAc_data_lit_review_v1.csv",
#                 na = "")

Check of table 1 from paper (selected sig proteins)

2025_12_24 - while most of the reanalysis matches the results presented in the paper. some did not quite. I suspected this was due to using the mouse database (which i excluded from my analysis). The code below is checking that.

# list of EE vs IC sig proteins shown in table 1 of the paper | 15 genes/proteins
acc_2_check <- c('P61751', 'Q64611',  'P24268', 'Q641Y8', 'Q63622', 'Q6PDL0',
                 'Q8BGY2', 'O54921', 'P28741', 'P21396', 'P36506', 'P10637',
                 'Q80ZA5', 'P07895', 'Q60930')

diff_check <- df_1 %>% 
  filter(Accession %in% acc_2_check)
head(diff_check, 13)
## # A tibble: 13 × 13
##    Accession Name        `Gene Symbol` `EC log2 fold` `EC est fold` `EC p value`
##    <chr>     <chr>       <chr>                  <dbl>         <dbl>        <dbl>
##  1 P61751    ARF4_RAT    ARF4                  -2.21          -4.63       0.0126
##  2 Q64611    CSAD_RAT    CSAD                  -5.30         -39.5        0.0203
##  3 P24268    CATD_RAT    CTSD                   1.32           2.49       0.0346
##  4 Q641Y8    DDX1_RAT    DDX1                  -0.622         -1.54       0.0203
##  5 Q641Y8    DDX1_RAT    DDX1                  -0.159         -1.12       0.758 
##  6 Q63622    DLG2_RAT    DLG2                   1.62           3.08       0.0328
##  7 Q63622    DLG2_RAT    DLG2                  -0.340         -1.27       0.990 
##  8 Q6PDL0    DC1L2_MOUSE DYNC1LI2               1.50           2.82       0.0429
##  9 Q8BGY2    IF5A2_MOUSE EIF5A2                 1.10           2.15       0.0446
## 10 Q8BGY2    IF5A2_MOUSE EIF5A2                 0.628          1.55       0.0874
## 11 O54921    EXOC2_RAT   EXOC2                  1.77           3.41       0.0328
## 12 P28741    KIF3A_MOUSE KIF3A                  2.23           4.68       0.0429
## 13 P21396    AOFA_RAT    MAOA                   0.818          1.76       0.0466
## # ℹ 7 more variables: `Coc log2 fold` <dbl>, `Coc est fold` <dbl>,
## #   `Coc p value` <dbl>, `interaction log2 effect size` <dbl>,
## #   `interaction effect size` <dbl>, `interaction p value` <dbl>,
## #   fraction <chr>
tail(diff_check, 13)
## # A tibble: 13 × 13
##    Accession Name        `Gene Symbol` `EC log2 fold` `EC est fold` `EC p value`
##    <chr>     <chr>       <chr>                  <dbl>         <dbl>        <dbl>
##  1 P21396    AOFA_RAT    MAOA                  0.818           1.76       0.0466
##  2 P21396    AOFA_RAT    MAOA                  0.0829          1.06       0.630 
##  3 P21396    AOFA_RAT    MAOA                  0.127           1.09       0.990 
##  4 P36506    MP2K2_RAT   MAP2K2               -2.25           -4.74       0.0181
##  5 P10637    TAU_MOUSE   MAPT                  1.22            2.33       0.0446
##  6 P10637    TAU_MOUSE   MAPT                  0.120           1.09       0.990 
##  7 Q80ZA5    S4A10_RAT   SLC4A10               3.02            8.10       0.0211
##  8 P07895    SODM_RAT    SOD2                  0.914           1.88       0.0487
##  9 P07895    SODM_RAT    SOD2                  0.851           1.80       0.114 
## 10 P07895    SODM_RAT    SOD2                  0.919           1.89       0.990 
## 11 Q60930    VDAC2_MOUSE VDAC2                 0.685           1.61       0.0300
## 12 Q60930    VDAC2_MOUSE VDAC2                 0.139           1.10       0.386 
## 13 Q60930    VDAC2_MOUSE VDAC2                 0.0991          1.07       0.990 
## # ℹ 7 more variables: `Coc log2 fold` <dbl>, `Coc est fold` <dbl>,
## #   `Coc p value` <dbl>, `interaction log2 effect size` <dbl>,
## #   `interaction effect size` <dbl>, `interaction p value` <dbl>,
## #   fraction <chr>

Bioinformatics resumed

2025_12_25/26 onwards

Running GO term GSEA with nested for loops

Currently only works on SBS Jennifer (insufficient memory on laptop?)

Previous Errors

  • Error in serialize(data, node$con) :… … error writing to connection
    • Fixed (kinda?) - Fine on SBS Jennifer server, issues is just on laptop
  • Error in `$<-.data.frame`(`*tmp*`, “.sign”, value = “activated”) : replacement has 1 row, data has 0
    • Fixed with if else for dotplots to only run if not empty (error was due to some GO results being empty due to no significant findings)
# not needed - works fine with uniprot entries
# # make trimmed entrez_id (rm everything after semi colon)
# df_4 <- df_3 %>%
#   ungroup() %>%
#   mutate(entrez_id = str_replace(
#     GeneID, pattern = ";.*", replacement = ""))

Setup containers & prepare info for sequences

# make vector of fractions (cytolsolic, membrane, & nuclear)
fraction_names <- unique(df_3$fraction)

# GO terms
GO_names <- c('BP', 'CC', 'MF', 'ALL')

# make container list to store results | empty for lists
gsea_GO_results_lists <- vector(mode = "list")

Operation

  1. For loop (for each fr fraction), compute a ranked genelist…

  2. Then nested for loop (for each go GO term),

    1. if, fr go result is not empty (> 0) then, make dot plot

    2. else, print message no result.

# works on SBS jennifer
for (fr in fraction_names){ # SEQUENCE - for each item (fr) in fraction_names
  
  # OPERATIONS
  # PART 1
  # make a ranked gene list for this fraction 
  geneList_working_df_go <- df_3 %>% # rats only mapped df from Lichti proteomics
    # filter by item | fraction is column in df_3, fr is i (current item in seq)

    filter(fraction == fr) %>%
    
    # compute ranking metric 
    mutate(ranking = sign(`EC log2 fold`)*(-log10(`EC p value`))) %>%
    # for any duplicate entries, # average & remove any nas
    # group_by(Gene.Names..primary.) %>%
    group_by(Entry) %>%
    summarise(ranking = mean(ranking, na.rm = TRUE)) %>%
    # lastly, arrange from largest positive to most negative
    arrange(-ranking)
  
  # Part 1.2
  # next, deframe into a named vector needed to supply to fgsea()
  geneList_vec_working_go <- tibble::deframe(geneList_working_df_go)
  
geneList_vec_working <- geneList_vec_working_go[!is.na(names(geneList_vec_working_go))]
  
  # store in list 
  gsea_GO_results_lists$geneList_vec[[fr]] <- geneList_vec_working_go
  
  # Part 1.3 | plots fine but does not save?
  # Plot ranked gene List
  gL_plot_go <- plot(geneList_vec_working_go, main = paste0(fr, '- Ranked gene list'))
  
  print(gL_plot_go) 
  
  gsea_GO_results_lists$geneList_plots[[fr]] <- gL_plot_go
  
  
  # PART 2 
  # nested for loop
  # so now for each fraction, go through all 3 go terms (BP, CC, MF)
  for (go in GO_names){ 
    
    # operation
    # run gseGO for each go term 
    GO_res_working <- clusterProfiler::gseGO(geneList = geneList_vec_working_go,
                                             OrgDb = org.Rn.eg.db,
                                             keyType = "UNIPROT",
                                             ont = go, # 
                                             minGSSize = 10,
                                             maxGSSize = 500,
                                             pvalueCutoff = 0.05,
                                             pAdjustMethod = "BH"
                                             )
    
  # save GO result original to list vec
  gsea_GO_results_lists$GO_results[[fr]][[go]] <- GO_res_working
    
    
  # extra (not really needed but save anyway
  gsea_GO_results_lists$GO_results_tib[[fr]][[go]] <- tibble::as_tibble(GO_res_working)
    

  # PART 3 - if else for dot plot
  if (nrow(GO_res_working@result) > 0 ) {
    
    #Make dot plots | this one saves properly
    # Part 3.1 - regular dot plots
    dp_go <- enrichplot::dotplot(GO_res_working,
                              showCategory = 30,
                              font.size = 11,
                              title = paste0(stringr::str_glue(
                                '{fr} - {go} - NAc - Lichti 2014 EE vs IC - gseGO'
                              )))
    print(dp_go) # print dot plots
    gsea_GO_results_lists$dotplots[[fr]][[go]] <- dp_go # save to list
    
  } else{print(paste0("No GO gsea result for dot plot in ", fr, "-", go, " :("))}
  
  } 
  
}

## NULL
## using 'fgsea' for GSEA analysis, please cite Korotkevich et al (2019).
## preparing geneSet collections...
## GSEA analysis...
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (33.82% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
## leading edge analysis...
## done...

## using 'fgsea' for GSEA analysis, please cite Korotkevich et al (2019).
## preparing geneSet collections...
## GSEA analysis...
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (33.82% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
## no term enriched under specific pvalueCutoff...
## [1] "No GO gsea result for dot plot in cytosolic-CC :("
## using 'fgsea' for GSEA analysis, please cite Korotkevich et al (2019).
## preparing geneSet collections...
## GSEA analysis...
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (33.82% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
## no term enriched under specific pvalueCutoff...
## [1] "No GO gsea result for dot plot in cytosolic-MF :("
## using 'fgsea' for GSEA analysis, please cite Korotkevich et al (2019).
## preparing geneSet collections...
## GSEA analysis...
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (33.82% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
## leading edge analysis...
## done...

## NULL
## using 'fgsea' for GSEA analysis, please cite Korotkevich et al (2019).
## preparing geneSet collections...
## GSEA analysis...
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (30.61% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize =
## minSize, : There were 22 pathways for which P-values were not calculated
## properly due to unbalanced (positive and negative) gene-level statistic values.
## For such pathways pval, padj, NES, log2err are set to NA. You can try to
## increase the value of the argument nPermSimple (for example set it nPermSimple
## = 10000)
## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize =
## minSize, : For some of the pathways the P-values were likely overestimated. For
## such pathways log2err is set to NA.
## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize =
## minSize, : For some pathways, in reality P-values are less than 1e-10. You can
## set the `eps` argument to zero for better estimation.
## leading edge analysis...
## done...

## using 'fgsea' for GSEA analysis, please cite Korotkevich et al (2019).
## preparing geneSet collections...
## GSEA analysis...
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (30.61% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize =
## minSize, : There were 7 pathways for which P-values were not calculated
## properly due to unbalanced (positive and negative) gene-level statistic values.
## For such pathways pval, padj, NES, log2err are set to NA. You can try to
## increase the value of the argument nPermSimple (for example set it nPermSimple
## = 10000)
## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize =
## minSize, : For some of the pathways the P-values were likely overestimated. For
## such pathways log2err is set to NA.
## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize =
## minSize, : For some pathways, in reality P-values are less than 1e-10. You can
## set the `eps` argument to zero for better estimation.
## leading edge analysis...
## done...

## using 'fgsea' for GSEA analysis, please cite Korotkevich et al (2019).
## preparing geneSet collections...
## GSEA analysis...
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (30.61% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize =
## minSize, : There were 1 pathways for which P-values were not calculated
## properly due to unbalanced (positive and negative) gene-level statistic values.
## For such pathways pval, padj, NES, log2err are set to NA. You can try to
## increase the value of the argument nPermSimple (for example set it nPermSimple
## = 10000)
## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize =
## minSize, : For some of the pathways the P-values were likely overestimated. For
## such pathways log2err is set to NA.
## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize =
## minSize, : For some pathways, in reality P-values are less than 1e-10. You can
## set the `eps` argument to zero for better estimation.
## leading edge analysis...
## done...

## using 'fgsea' for GSEA analysis, please cite Korotkevich et al (2019).
## preparing geneSet collections...
## GSEA analysis...
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (30.61% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize =
## minSize, : There were 25 pathways for which P-values were not calculated
## properly due to unbalanced (positive and negative) gene-level statistic values.
## For such pathways pval, padj, NES, log2err are set to NA. You can try to
## increase the value of the argument nPermSimple (for example set it nPermSimple
## = 10000)
## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize =
## minSize, : For some of the pathways the P-values were likely overestimated. For
## such pathways log2err is set to NA.
## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize =
## minSize, : For some pathways, in reality P-values are less than 1e-10. You can
## set the `eps` argument to zero for better estimation.
## leading edge analysis...
## done...

## NULL
## using 'fgsea' for GSEA analysis, please cite Korotkevich et al (2019).
## preparing geneSet collections...
## GSEA analysis...
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (89.71% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
## no term enriched under specific pvalueCutoff...
## [1] "No GO gsea result for dot plot in membrane-BP :("
## using 'fgsea' for GSEA analysis, please cite Korotkevich et al (2019).
## preparing geneSet collections...
## GSEA analysis...
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (89.71% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
## no term enriched under specific pvalueCutoff...
## [1] "No GO gsea result for dot plot in membrane-CC :("
## using 'fgsea' for GSEA analysis, please cite Korotkevich et al (2019).
## preparing geneSet collections...
## GSEA analysis...
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (89.71% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
## no term enriched under specific pvalueCutoff...
## [1] "No GO gsea result for dot plot in membrane-MF :("
## using 'fgsea' for GSEA analysis, please cite Korotkevich et al (2019).
## preparing geneSet collections...
## GSEA analysis...
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (89.71% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
## no term enriched under specific pvalueCutoff...
## [1] "No GO gsea result for dot plot in membrane-ALL :("
# # horrific name but at least clear...
# Lichti_2014_gsea_GO_results_2025_12_26_v1 <- gsea_GO_results_lists
# 
# save(Lichti_2014_gsea_GO_results_2025_12_26_v1,
#      file = "output/2025_12_26_Lichti_2014_gsea_GO_results_rats_NAc_v1.rda")
# cow plot for combined + plotlist = to plot a list of plots 
# cowplot::plot_grid(plotlist = c(gsea_GO_results_lists$dotplots),
#                    ncol = 2)

# cow plot not working 
# Warning message:
# In as_grob.default(plot) : Cannot convert object of class list into a grob.


# # this works 
# print(gsea_GO_results_lists$dotplots)
sessionInfo()
## R version 4.5.2 (2025-10-31)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 22.04.5 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0  LAPACK version 3.10.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Pacific/Auckland
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] enrichplot_1.28.4      clusterProfiler_4.16.0 org.Rn.eg.db_3.21.0   
##  [4] AnnotationDbi_1.70.0   IRanges_2.42.0         S4Vectors_0.46.0      
##  [7] Biobase_2.68.0         BiocGenerics_0.54.0    generics_0.1.4        
## [10] UniProt.ws_2.48.0      gt_1.2.0               lubridate_1.9.4       
## [13] forcats_1.0.1          stringr_1.5.2          dplyr_1.1.4           
## [16] purrr_1.1.0            readr_2.1.5            tidyr_1.3.1           
## [19] tibble_3.3.0           ggplot2_4.0.0          tidyverse_2.0.0       
## [22] readxl_1.4.5          
## 
## loaded via a namespace (and not attached):
##   [1] DBI_1.2.3               gson_0.1.0              httr2_1.2.1            
##   [4] AnVILBase_1.2.0         rlang_1.1.6             magrittr_2.0.4         
##   [7] DOSE_4.2.0              compiler_4.5.2          RSQLite_2.4.5          
##  [10] png_0.1-8               vctrs_0.6.5             reshape2_1.4.4         
##  [13] pkgconfig_2.0.3         crayon_1.5.3            fastmap_1.2.0          
##  [16] dbplyr_2.5.1            XVector_0.48.0          labeling_0.4.3         
##  [19] utf8_1.2.6              rmarkdown_2.30          tzdb_0.5.0             
##  [22] UCSC.utils_1.4.0        bit_4.6.0               xfun_0.53              
##  [25] cachem_1.1.0            aplot_0.2.9             GenomeInfoDb_1.44.3    
##  [28] jsonlite_2.0.0          progress_1.2.3          blob_1.2.4             
##  [31] BiocParallel_1.42.2     parallel_4.5.2          prettyunits_1.2.0      
##  [34] R6_2.6.1                bslib_0.9.0             stringi_1.8.7          
##  [37] RColorBrewer_1.1-3      jquerylib_0.1.4         cellranger_1.1.0       
##  [40] GOSemSim_2.34.0         Rcpp_1.1.0              knitr_1.50             
##  [43] ggtangle_0.0.9          R.utils_2.13.0          BiocBaseUtils_1.10.0   
##  [46] igraph_2.1.4            splines_4.5.2           Matrix_1.7-4           
##  [49] timechange_0.3.0        tidyselect_1.2.1        qvalue_2.40.0          
##  [52] rstudioapi_0.17.1       yaml_2.3.10             codetools_0.2-20       
##  [55] curl_7.0.0              rjsoncons_1.3.2         lattice_0.22-7         
##  [58] plyr_1.8.9              treeio_1.32.0           withr_3.0.2            
##  [61] KEGGREST_1.48.1         S7_0.2.1                evaluate_1.0.5         
##  [64] gridGraphics_0.5-1      BiocFileCache_2.16.2    xml2_1.4.0             
##  [67] Biostrings_2.76.0       ggtree_3.17.0           pillar_1.11.1          
##  [70] filelock_1.0.3          ggfun_0.2.0             hms_1.1.3              
##  [73] tidytree_0.4.6          scales_1.4.0            glue_1.8.0             
##  [76] lazyeval_0.2.2          tools_4.5.2             data.table_1.17.8      
##  [79] fgsea_1.34.2            fs_1.6.6                fastmatch_1.1-6        
##  [82] cowplot_1.2.0           grid_4.5.2              ape_5.8-1              
##  [85] nlme_3.1-168            patchwork_1.3.2         GenomeInfoDbData_1.2.14
##  [88] cli_3.6.5               rappdirs_0.3.3          gtable_0.3.6           
##  [91] R.methodsS3_1.8.2       yulab.utils_0.2.1       sass_0.4.10            
##  [94] digest_0.6.37           ggrepel_0.9.6           ggplotify_0.1.3        
##  [97] farver_2.1.2            memoise_2.0.1           htmltools_0.5.8.1      
## [100] R.oo_1.27.1             lifecycle_1.0.4         httr_1.4.7             
## [103] GO.db_3.21.0            bit64_4.6.0-1